#clean up
rm(list=ls())

#Set Working Directory
#setwd('/Volumes/MONOGAN/stateAreal/')
#setwd("C:/Users/Josh Jackson/Dropbox/Spatial Paper with Jamie/Replication")

####Decision Tree Graph, Figure 1####
#postscript("decisionTree.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("decisionTree.pdf",width=5, height=5, pointsize=12)
pdf("figure1.pdf",width=5, height=5, pointsize=12)
plot(x=c(0,1.25),y=c(0,1),xlab="",ylab="",axes=F,type="n")
arrows(x0=.5,y0=.9,x1=.25,y1=.5,length=.1)
arrows(x0=.5,y0=.9,x1=.75,y1=.5,length=.1)
text(x=.5,y=.9,label="Functional Form",pos=3)
text(x=.35,y=.75,label="Diffusive",pos=2)
text(x=.65,y=.75,label="Confined",pos=4)
text(x=.25,y=.5,label="Spatial Lag Model",pos=1)
text(x=.775,y=.5,label="Error Structure",pos=1)
arrows(x0=.75,y0=.4,x1=.6,y1=.1,length=.1)
arrows(x0=.75,y0=.4,x1=.9,y1=.1,length=.1)
text(x=.675,y=.25,label="Direct",pos=2)
text(x=.825,y=.25,label="Hierarchical",pos=4)
text(x=.5,y=.1,label="MLE CAR/SAR",pos=1)
text(x=1,y=.1,label="Bayesian CAR",pos=1)
dev.off()


####Forest Plot Results####
#Cross-Sectional Results, draw Figure 4
#EWM RESULTS:
#postscript("ewmForest.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("ewmForest.pdf",width=5, height=5,  pointsize=10)
pdf("figure4.pdf",width=5, height=5,  pointsize=10)
lin.mean<-c(.7626, .8414, -1.023, .1294, -.2245, .7866, 1.053, .6601, .4428, -.3308, .4775)
lin.lower<-c(.5964, .5443, -1.32, -.1139, -.4867, .6171, .97, .5772, .2218, -.4606, .2508)
lin.upper<-c(.9294, 1.138, -.7259, .372, .03814, .9564, 1.136, .7431, .6647, -.2013, .7049)

car.mean<-c(0.6033, 0.8221,-0.8485,0.06182,-0.1247,0.6976,1.066,0.6652,0.3687,-0.2959,0.4933)
car.lower<-c(0.4052,0.4937,-1.215,-0.2279,-0.4378,0.5007,0.9106,0.5215,0.09225,-0.4948,0.2178)
car.upper<-c(0.8015,1.152,-0.4791,0.3506,0.1898,0.8948,1.221,0.8088,0.6458,-0.0954,0.7687)

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
   #CAR credible intervals wider across the board

par(omi=c(.2,1.1,.1,.01))
plot(x=lin.mean,y=c(1:11-.1),xlim=c(-1.3,1.3),ylim=c(0,12), xlab="Coefficient", ylab="",axes=F)
axis(1)
axis(2, at=c(1:11), labels=c('Opinion on party elite','Opinion on Dem ID','Party elite on Dem ID','Opinion on strength','Party elite on strength','Dem ID on strength','Party elite on legis lib','Strength on legis lib','Opinion on policy','Strength on policy','Legis lib on policy'),las=1)
box()
abline(v=0,col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(1:11)-.1,y1=c(1:11)-.1, lty=2)
points(x=car.mean,y=c(1:11+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:11+.1),y1=c(1:11+.1), col='red')
legend(x=-1.3, y=2.5, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))
dev.off()

#Duration Results, draw Appendix Figure A.4
#Lottery RESULTS:
#postscript("lotteryForest.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("lotteryForest.pdf",width=5, height=5,  pointsize=10)
pdf("figureA4.pdf",width=5, height=5,  pointsize=10)
lin.mean<-c(.6096,.504,-.004703,.1593,-.1271,.5193,.05508,-.02823)
lin.lower<-c(-.1845,-.224,-.1197,-.4473,-.1786,.2399,.008679,-.04531)
lin.upper<-c(1.43,1.279,.1307,.7681,-.0763,.7999,.09814,-.0120)

car.mean<-c(0.5633,.5017,.05112,0.1383,-.103,0.08812,.06677,-.03728)
car.lower<-c(-.2246,-.2354,-.0601,-.4904,-.1584,-.2188,.0249,-.05607)
car.upper<-c(1.388,1.28,.1593,0.7553,-.05092,0.3978,.1095,-.01963)

par(omi=c(.2,.9,.1,.01))
plot(x=lin.mean,y=c(1:8-.1),xlim=c(-.75,1.5),ylim=c(0,11), xlab="Coefficient", ylab="",axes=F)
axis(1)
axis(2, at=c(1:8,10), labels=c('Election 1','Election 2','Lagged income','Unified party','Fundamentalists','Neighbors adopted','Time','Time squared', 'Lagged fiscal health'),las=1)
box()
#abline(v=0,col='gray60')
segments(x0=0,x1=0,y0=-1,y1=8.6, col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(1:8-.1),y1=c(1:8-.1), lty=2)
points(x=car.mean,y=c(1:8+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:8+.1),y1=c(1:8+.1), col='red')
#legend(x=-1, y=1, pch=c(2,1), lty=c(2,1), col=c('red','black'), legend=c('CAR', 'Independent'))
abline(h=8.6,col='black')
abline(h=8.7,col='black')

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
#mean(abs(car.upper-car.lower)-abs(lin.upper-lin.lower))
   #CAR wider on 5 out of 8. 

par(new=T,omi=c(.2,.9,.1,.01))
lin.mean<-c(-3.422)
lin.lower<-c(-6.847)
lin.upper<-c(-.2447)

car.mean<-c(-3.343)
car.lower<-c(-6.817)
car.upper<-c(-.08333)

abs(car.upper-car.lower)-abs(lin.upper-lin.lower)
   #inconsistent
	#CAR wider on 6 out of 9.

plot(x=lin.mean,y=c(10-.1),xlim=c(-9,4),ylim=c(0,11), xlab="Coefficient", ylab="",axes=F)
axis(3)
segments(x0=0,x1=0,y0=8.7,y1=12, col='gray60')
segments(x0=lin.lower,x1=lin.upper,y0=c(10-.1),y1=c(10-.1), lty=2)
points(x=car.mean,y=c(10+.1), col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(10+.1),y1=c(10+.1), col='red')
legend(x=-1.75, y=8.5, pch=c(2,1), lty=c(1,2), col=c('red','black'), legend=c('CAR', 'Independent'))
dev.off()


#Multilevel model results (Margalit 2013), draw Figure 5 
#postscript("margalitForest.eps", width=5, height=5, paper='special', onefile=F, pointsize=10)
#pdf("margalitForest.pdf", width=5, height=5, pointsize=10)
pdf("figure5.pdf", width=5, height=5, pointsize=10)
sw.mean=   c(.041,-.006, .000, .018, .171,-.122, .003,-.005, .009,-.008)
sw.lower=  c(.024,-.022,-.017, .003, .146,-.145,-.013,-.020,-.021,-.025)  
sw.upper=  c(.060, .010, .017, .033, .195,-.098, .019, .009, .038, .008)

car.mean=  c(.042,-.006, .000, .018, .171,-.122, .003,-.005, .009,-.006)
car.lower= c(.024,-.022,-.018, .003, .146,-.145,-.013,-.020,-.021,-.021)
car.upper= c(.060, .010, .018, .033, .195,-.098, .019, .010, .038, .010)

par(omi=c(.1,.8,.1,.1))
plot(x=sw.mean, y=c(1:10)-.15, ylim=c(0,11.5), xlim=c(-.15,.2), xlab="Coefficient", ylab="", axes=F)
axis(1,cex.axis=.8)
axis(2, at=c(1:11), labels=c("Lost Job", "HH Income Drop", "Job Less Secure", "Sp. Lost Job", "Democrat", "Republican", "Long Term Unemp.", "New Job", "Not Looking", "Unemployment", "Prev. Welfare"), las=1)
box()
#abline(v=0, col="gray60")
segments(x0=sw.lower,x1=sw.upper,y0=c(1:10)-.15, y1=c(1:10)-.15, lty=2)
points(x=car.mean,y=c(1:10)+.15, col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:10)+.15,y1=c(1:10)+.15, col='red')
abline(h=10.6, col="black")
abline(h=10.5, col="black")
segments(x0=0,x1=0,y0=-1,y1=10.5, col='gray60')
legend(x=0.05, y=9, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))

#"Prev. Welfare"
par(new=T, omi=c(.1,.8,.1,.1))
sw.mean2= .588
sw.lower2= .569
sw.upper2= .608

car.mean2= .589
car.lower2= .569
car.upper2= .608

plot(x=sw.mean2,y=c(11-.15),ylim=c(0,11.5),xlim=c(0,0.65), xlab="Coefficient", ylab="",axes=F)
axis(3, cex.axis=.8)
segments(x0=0,x1=0,y0=10.6,y1=12, col='gray60')
segments(x0=sw.lower2,x1=sw.upper2,y0=c(11-.15),y1=c(11-.15), lty=2)
points(x=car.mean2,y=c(11+.15), col='red',pch=2)
segments(x0=car.lower2,x1=car.upper2,y0=c(11+.15),y1=c(11+.15), col='red')
dev.off()

abs(car.upper-car.lower)-abs(sw.upper-sw.lower)
   #The CAR interval is wider in two cases, and the independent is wider in one. In all other cases, the intervals are equivalent to our level of precision.
abs(car.upper2-car.lower2)-abs(sw.upper2-sw.lower2)
   #Identical 

#Clean Air results, draw Appendix Figure A.6
#postscript("cleanairForest.eps", width=5, height=5, paper='special', onefile=F, pointsize=10, family='ComputerModern')
#pdf("cleanairForest.pdf", width=5, height=5, pointsize=10)
pdf("figureA6.pdf", width=5, height=5, pointsize=10)

ind.mean= c(.0709, .0572, .0222, .0040, .1836, .0278, .0014)
ind.lower=c(.0624, .0502, .0145,-.0023, .1667, .0230,-.0046)
ind.upper=c(.0792, .0640, .0299, .0101, .1999, .0325, .0073)

car.mean= c(.0705, .0573, .0219, .0038, .1858, .0278, .0015)
car.lower=c(.0604, .0488, .0127,-.0037, .1642, .0219,-.0058)
car.upper=c(.0804, .0657, .0313, .0112, .2069, .0335, .0089)

par(omi=c(.1,.8,.1,.1))
plot(x=ind.mean, y=c(1:7)-.15, ylim=c(0,7.5), xlim=c(-.015, .215), xlab="Coefficient", ylab="", axes=F)
axis(1)
axis(2, at=c(1:7), labels=c("Fed. Actions (lag)", "Dems Unified", "Reps Unified", "Unemployment", "Manufacturing", "Grants", "Groups"), las=1)
box()
segments(x0=ind.lower,x1=ind.upper,y0=c(1:7)-.15, y1=c(1:7)-.15, lty=2)
points(x=car.mean,y=c(1:7)+.15, col='red',pch=2)
segments(x0=car.lower,x1=car.upper,y0=c(1:7)+.15,y1=c(1:7)+.15, col='red')
abline(v=0, col='gray60')
legend(x=0.12, y=1, pch=c(2,1), lty=c(1,2), col=c('red', 'black'), legend=c('CAR', 'Independent'))
dev.off()

abs(car.upper-car.lower)-abs(ind.upper-ind.lower)
#wider across the board

####RAW DATA PLOTS####
#Lottery Data, draw Appendix Figure A.3
#setwd("/Volumes/MONOGAN/stateAreal/lottery/")
bbNew<-read.csv("bbExtended4.csv")
bbNew$zInc<-(bbNew$lagIncome-mean(bbNew$lagIncome))/sd(bbNew$lagIncome)
bbNew$zFisc<-(bbNew$lagFiscal-mean(bbNew$lagFiscal))/sd(bbNew$lagFiscal)
bbNew$zRelig<-(bbNew$religion-mean(bbNew$religion))/sd(bbNew$religion)
lotto<-bbNew
table(lotto$year)
count<-abs(table(lotto$year)-50)[-1]
count<-c(0,count,43,43,43,43,44)
years<-c(1963:2017)
names(count)<-years

#postscript("lotto.eps",horizontal=FALSE, width=5, height=5, paper='special', onefile=FALSE, pointsize=10)
#pdf("lotto.pdf",width=5, height=5, pointsize=10)
pdf("figureA3.pdf",width=5, height=5, pointsize=10)
plot(y=count,x=years,xlab="Year",ylab="Number of States With Lotteries",type='o',ylim=c(0,50),lwd=2,pch=20)
abline(h=0,col='gray60')
dev.off()


#Clean Air Enforcement Data, draw Appendix Figure A.5
library(lattice)
library(foreign)
#setwd("/Volumes/MONOGAN/stateAreal/wood1992/")
caa<-read.dta('cleanAir0109.dta')
caa<-subset(caa, year>=2001)
caa<-subset(caa, year<=2009)

#trellis.device("pdf",color=FALSE,file="stateActions0109.pdf")
#trellis.device("pdf",color=TRUE,file="stateActions0109.pdf")
#postscript("stateActions.eps", horizontal=FALSE, colormodel='rgb', width=5, height=5, paper="special", onefile=FALSE, family='ComputerModern')
trellis.device("pdf",color=FALSE,file="figureA5.pdf")
xyplot(log(stateact)~year, caa, type='l', groups=state, xlab="Year", ylab="State Actions (Logged Scale)")
dev.off()
